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Abstract 

Energy loss through optically thin radiative cooling plays an important part 
in the evolution of astrophysical gas dynamics and should therefore be consid- 
ered a necessary element in any numerical simulation. Although the addition 
of this physical process to the equations of hydrodynamics is straightforward, 
it does create numerical challenges that have to be overcome in order to en- 
sure the physical correctness of the simulation. First, the cooling has to 
be treated (semi-)implicitly, owing to the discrepancies between the cooling 
timescale and the typical timesteps of the simulation. Secondly, because of 
its dependence on a tabulated cooling curve, the introduction of radiative 
cooling creates the necessity for an interpolation scheme. In particular, we 
will argue that the addition of radiative cooling to a numerical simulation 
creates the need for extremely high resolution, which can only be fully met 
through the use of adaptive mesh refinement. 
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1. Introduction 

A gas at high temperature loses energy through radiation. In astrophysics 
the most common process to produce this radiation is through the recombi- 
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nation of ionized particles. This process can cause a significant decrease in 
local temperature over a short period of time, depending on the local con- 
ditions. In order to produce a correct numerical simulation of the evolution 
of a gas-structure on astrophysical scale, a method to quantify this form of 
radiative cooling has to be implemented. The most simple form is that of 
optically thin radiative cooling. Here it is assumed that the gas in which the 
photons are emitted is completely optically thin, so that any photon that is 
emitted will simply leave the physical domain of the simulation rather than 
being absorbed elsewhere. This is a valid approach for astrophysical phenom- 
ena, which tend to have low to extremely low densities. Indeed, the typical 
particle density of the interstellar medium is on the order of io~ 24 - -23 g/cm 3 
(1...10 particles per cm 3 ), whereas the gas density of earth's atmosphere at 
sea level is approximately 1.2 x 10 -3 g/cm 3 . The great advantage of the 
optically thin approach over a more complicated radiative transfer method is 
that it reduces the physics of radiative cooling to a purely local phenomenon. 
The energy of the gas in a given spot decreases, owing to its local density, 
temperature and degree of ionization. This maintains the localized nature of 
the Euler equations of gas dynamics (eq. [TJ. Also, since the energy content of 
the escaping photons is lost to the gas anyway, there is no need to introduce 
an extra conserved variable for the radiation field internal energy. 

The behaviour of an ideal gas can be described by solving the conse rvation 



equations for mass, momentum and energy here in Euleria n form (e.g. lLaney 
19981 : ICastorl . Il998l ; iGoedbloed. Keppens fc Poedtsl . l2010l ) 

dp 



at + v ' ( " v) = 



<9(pv) 



()l V-(pvv) + Vp = F, (1) 

— + V • (ev + pv) = G + v • F, 
at 

with p the mass density, v the velocity vector, e the energy density, p the 
thermal pressure and F and G the momentum and energy source terms. 
These source terms can generally include such effects as radiation pressure, 
gravity etc. In our simulations, we will use optically thin radiative cooling 
as the only source term, making F = and G as specified in eq. HI Since the 
pressure and density are related through 

e = ^ + J^Ty 0) 
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with 7 the adiabatic constant, this set of equations forms a closed system 
that can be solved for a set of boundar y values. 



O ptically thin radia tive cooling (e.g. lMacDonald fc Baileyi . ll98ll ; lMeUema &: Lundqvist 



20021 ; iTownsendl l2009h is given by the equation 
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with dp/dt the change in thermal pressure due to the cooling, rij and n e 
are the local ion particle density and electron density respectively, which 
relate to the total gas density through p = riirrii + n e m e , with m e and 
TOj the electron and ion masses. A(T) is a cooling function dependent 
on the local temperature T and is usually taken from a table that is ei- 
ther constructed through observation, combining insights from astrophysics 
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2002|; 
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2009 


. In astrophysical literature 



one can find different tables, being constructed and gradually improved upon, 
that depend on the composition of the gas. The function C in the alternative 
way of writing the source term in eq. |3] specifies the energy losses per unit 
mass. 

1.1. Numerical challenges 

The cooling is implemented in a gas dynamics code by adding the effect 
of the right-hand term in eq. |3] to the Euler equation for conservation of 
total energy (eq. [1]), which then includes an energy sink term G and takes 
the form 

de 

— + V-(ev) + V-(pv) = -mn e A(T). (4) 

Though this is simple enough by itself, the unique nature of the radiative 
cooling function may lead to some difficulty with the implementation. Specif- 
ically, depending on the local characteristics of the gas, the cooling timescale 



7~cool 
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n;n e A(T) 



(5) 



can be much shorter than the dynamical timescale of the gas itself. Therefore, 
additional modifications to the original code are required in order to either 
limit the timestep to the cooling timescale (Eq. |5]), or use an implicit calcu- 
lation scheme for the cooling, even if the gasdynamics themselves are solved 
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explicitly. The latte r method, which has the advantage of st ability inherent 
in implicit methods ( iBuelow. Venkateswaran fc Merkld . l200ll ). is usually pre- 
ferred as the cooling timescales can become very short indeed, which would 
slow down numerical simulations to an unacceptable level. Moreover, the 
typically tabulated nature of the cooling function A(T) introduces the need 
for adequate interpolation formulae as well as potential difficulties across 
temperature ranges with strong variations (since the derivative A'(T) is not 
known analytically). 

1.2. Astrophysical challenges 

Taking radiative cooling into account can change the result of a numerical 
simulation to a considerable extent. First of all, depending on circumstances, 
a large part of the available energy can 'leak' out of the gas through radiation. 
If this happens it can change the morphology of the g of high 

thermal pressure disappear and the gas makes a transition from adiabatic to 
isothermal behavior. 

A second effect is the sudden formation of radiative cooling instabilities. 
These can result from either the density or the temperature dependence of the 
cooling. In the case of density dependent cooling instabilities, high density 
regions lose more energy than their surroundings, and therefore have a lower 
temperature, which leads to a lower local pressure. Hence, these regions will 
be compressed, causing an increase in local density, which in turn means 
a further increase in the cooling rate. This process can repeat itself until 
the temperature of the gas drops to zero, or whatever floor temperature has 
been set for the simulation, depending on the physical conditions one tries to 
simulate. As a result, simulations that include radiative cooling tend to show 
high density regions of much smaller physical size and higher local density 
than if the same model had been run without radiative cooling. 

Thermal instabilities can also result under isobaric conditions, in those 
temperature ranges where A(T) decreases for increasing T. A runaway con- 
densation results when 
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with C the energy gain/loss per unit mass due to non- adiabatic processes. 
Similarly, sound waves can turn into overstable m odes under isentropic con- 
ditions [8C/dT]s < fjFieldl . Il951t iParkerl . Il953l ). where S = pp" 1 denotes 
the entropy. 
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Table 1: Physical input parameters 

10" b M /yr (= 6.3 x 10 19 g/s) 
Ko 1500 km/s 

Pism 10~ 22 - 5 g/cm 3 

Tism 100 K 



Moreover, the high density structures formed by the cooling can also 
serve as seed for other instabilitie s, such as Kely i n-Hel mhortz (when two 



parallel flows have a strong shear (1 Chandr asekharl . Il98ll ) ) , Rayleigh- Taylor 
(whe n gravitational, c entrifugal or thermal pressure driven accel eration oc- 
curs Il98ll : iFournier. Gauthier fc Renaudl . 120021 )) and thin 
shell instabilities (which are the result of a thin shell being compress ed be- 
tween two areas with asymmetric pressure gradients (IVishniad . Il983l a As 
we will demonstrate in this paper, resolving such high density structures 
presents a challenge to the numerical code and necessitates an increase in 
resolution, which can best be achieved through the use of adaptive mesh 
refinement. 

Finally, radiative cooling allows us to compare our simulations directly 
with observations, as the radiative flux produced in this way is an observable 
quantity. 



2. Astrophysical application: stellar wind expansion 

As our test case, we use the expansion of a supersonic stellar wind into 
a surrounding constant densi ty interstellar medium (ISM ) . Thi s problem was 
appro xi mated analytically bvlCastor. McCray fc Weaver! (119751 ) , IWeaver et al 



( 119771 ) , lOstriker fc McKed (119881) and others and has since then been tested 



numerically by e.g. iGarcfa-Segura. Mac Low &: Langerl ( 1l996al lbl). with the 
numerical results showing close qualitative and quantitative agreement with 
the original analytical model. 

Schematically, the wind blown bubble shows the following structure. The 
stellar wind collides with the ISM. As the wind slows down because of the 
collision, its kinetic energy is converted to thermal energy, creating a 'hot 
bubble' of shocked wind material, which is contained by the reverse shock 
(Rl) on one end and a contact discontinuity (R2) on the other. At the contact 
discontinuity (R2), the high thermal pressure of the 'hot bubble' allows it 
to expand outward into the ISM, sweeping up a shell of shocked ISM (see 
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fig. H] and further). According to the analytical model the outer radius of the 
shocked ISM shell should be at time t at position: 



R?, 



250 
308tt 



1/5 



r 1/5 -1/5,3/5 



(7) 



with L w = \^jfV£ the mechanical luminosity of the wind with dM/dt the 
mass loss rate, pism the density of the ambient i nterstellar me d ium, the 



19770 . The 



terminal velocity of the wind and t the time ([Weaver et al 
numerical factor is a result of the choice of units, in this case cgs (centimeter, 
gram, second). 

It should be noted here that this model is somewhat simplified from an 
astrophysical point of view, as it assu mes that the stellar wind expands into 
a cold medium. More recent models (IFreyer. Hensler. fc Yorkd . 120031 12006 ; 
van Marie. Langer fc Garcia-Segural . 12005k 120071 ) , which include the effect of 
photo-ionization by radiation from the progenitor star show a rather more 
complicated result. In fact, even these models neglect external factors such 
as stellar winds and ionizing radiation from neighboring star systems. Such 
facto rs cause inhomogeneit ies in the ISM, which complicate the matter fur- 
ther (IMellema et al.l . 120061 ). Nevertheless, this is a good test case, since it 
is extremely well documented and the result can change considerably due to 
the influence of radiative cooling. For our input parameters we choose the 
parameters shown in Table 1. The mechanical luminosity (L w ) which follows 
from these parameters is high, with the mass loss rate dMj dt and wind ve- 
locity reflecting values typically found only in extr emely massiv e stars 
(> 60 71^) of O and hydrogen rich Wolf-Rayet type (ILanger et al.l. Il994l: 
Garcia-Segura. Mac Low fc Langerl . Il996al ; Ivan Marie. Langer fc Garcia-Segura , 
20071 ). These values were chosen deliberately to create a powerful, high den- 
sity shock which causes a high cooling rate. 



3. Numerical approach 



For our simulations we use the AMRVAC code ( iKeppens et al.l . 12003 



van der Hoist fc Keppend . 120071 ). This is a fully conservative code, which 



can solve e.g. the Euler equations of hydrodynamics on a grid that can 
be either Cartesian, cylindrical or spherical. The code includes an adaptive 
mesh refinement (AMR) scheme to adjust the resolution depending on preset 
criteria. In addition, AMRVAC allows the use of a wide variety of shock 
capturing solvers. 
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Figure 1: Cooling table for gas at solar metallicity according to iMellema fc Lundqvist 
(2002). Note the steep jump in cooling rate around 10 4 K, which corresponds to the 
(de-)ionization of hydrogen. 
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3.1. Discretization of the Euler equations 

The numerical form of the E uler equations, as used in the AMRVAC code 



has been described before (e.g. iKeppens et al.l . 120031 ). but we briefly recall 



the main solver strategy as employed for the simulations reported here. The 
integration scheme advances conserved variables (mass, momentum and total 
energy), which represent grid cell values in the standard finite volume sense. 
The numerical update takes the general form (specified to ID Cartesian for 
simplicity) 



predictor 
corrector 



u, 



n+1 



u? - 



ui 



2Ax ( F[U ^> 



~F(U^ /2 , 



A* 
Ac 



r i+l/2 



LF,n+| 
r i-l/2 I ■ 



with Ui the set of conserved variables (p, pv,pv 2 /2 +p/( i y — 1)) at the cell 
center of grid i, and F the physical fluxes deduced from eq. [T]to be (pv, pv 2 + 
p, ev +pv), which are in the (non-conservative) Hancock predictor step eval- 



uated at cell interface states U 



L,n 



u 



R,n 



i+i/2i ■'■ n ^ ne corrector step, we use the 

conservative To tal Variation Diminishing Lax-Friedrichs (TVDLF) (jToth &: Odstrcil 
19961 ; lYed . 119891 1 spatial discretization. This method specifies the fluxes F^, 2 
according to 



r i+l/2 



F(U t L + m) + F(U* 1/2 ) 



+1/2) ' 1 \ u i+y 

■L _j_ TjR 
max/ i+l/2^ U i+l/2 



TjR 

U i+l/2 



U, 



i+1/2 



(9) 



The locally computed c max denotes the maximum physical propagation speed 
at the (averaged) cell interface state, and is for hydrodynamics given by 
M + \fyp/p- To obtain the cell interface states from the cell center values, 
a limited linear reconstruction determines the Uk_ 1/2 and UR l/2 states as 



U i+l/2 

TjR 

U i+l/2 



+ 1/2 

Ut + AUt/2, 
U i+1 - Af7 m /2, 



(10) 



which involves a limited slope At/, = AUi<j)(ri) = AUi<j)(AUi-i / AUi) . The 
limiter is here written to act on the cell differences AU = U + \ — Ui, but 
in AMRVAC these can also be employed on the corresponding primitive 
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variables (p,v,p) or other user selected (e.g. logarithmically stretched) com- 
binations. Many flavors are implemented, but for most of our calculations 
we use the 'minmod' flux limiter, which gives for </>(r) the behavior 

1 if r > 1, 

r ifO>r>l, (11) 
if r < 0. 

To compare the effect of differe nt flux hunters on the result, we also run 



simulations with the 'van Leer' ( Ivan Leerl . 11974 ) limiter 



ix- if>>0, 
if r < 0. 

Since our astrophysical problem involves density and pressure jumps of sev- 
eral orders of magnitude, we use the logarithm of these positive primitive 
variables in the limiter. 

In order to avoid numerical instability in this explicit two-step integra- 
tion procedure, the size of each timestep has to be limited by the Courant- 
Friedrichs-Levy (CFL) condition, which stipulates that the time step must 
obey 

Ax 



At n < 



(13) 



Here Ax is the (local when AMR) size of a gridcell. The current AMRVAC 
code uses a single timestep for the entire grid hierarchy in the AMR, so the 
righthand value of Eq. [T3] is calculated for all gridcells to find the global 
minimal value. The new timestep becomes a fixed fraction of this value, in 
our case 0.25. In the simulations discussed below, this scheme is only slightly 
modified to handle the actual spherical geometry, and multi-dimensional sim- 
ulations handle fluxes from multiple directions in unsplit fashion. 

3.2. Radiative cooling methods 

In order to simulate the effect of radiative cooling, an extra module has 
been added to the code, which updates the local energy of the gas according 
to Eq. [3j This new module allows for various approaches to handle the sink 
term. The simplest method is a fully explicit cooling routine: 

e n+l = e n _ n ^A(T n )At n , (14) 
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where n and n + 1 denote discrete time levels before and after the timestep 
At n respectively. For this approach, numerical stability is maintained by 
forcing the timestep At n to remain below a pre-set fraction of the cooling 
timescale (Eq. [5]) everywhere in the grid, this upper limit being enforced on 
top of the standard CFL condition. 

Alternatively, the cooling can be calculated using a semi-implicit scheme 
so that: 

e n+l = e n _ n « n ™ A (T n+1 )At n . (15) 

This way, no extra limit has to be placed on the size of the timestep, but 
computationally it is expensive, since the numerical scheme must estimate 
the cooling rate at T n+l . In our implementation the cooling value A(T n+1 ) is 
estimated through a half-step refinement routine that iterates until a preset 
precision is achieved. This is a reliable method but the number of itera- 
tions can become large (typically about 10. ..20 iterations to reach a relative 
precision of 10~ 5 ). 

Several hybrid solutions have also been implemented. One uses a cooling 
value found by taking the average between A(T n ) and A(T n+1 ), without 
adjusting the size of the timestep, while the second allows the timestep to 
be split up into smaller, explicitly treated 'cooling steps'. The first has 
the advantage of providing a very fast solution at the expense of physical 
accuracy, whereas the latter is more precise, but can be nearly as slow as a 
fully explicit method, depending on the nature of the simulation. 

Finally, we have implemented the new exact-integration method proposed 



by iTownsendl (120091 ) . This method uses an exact integration of the cooling 
equation [3j The temperature evolution due to radiative cooling of a gas 
parcel at constant density is calculated in advance, starting at extremely 
high temperature (usually 10 8 K) and followed all the way down to the lowest 
temperature (usually 10 2 K) for which a cooling value is known. Since this is 
only done once and does not have to be repeated during the actual simulation, 
it can be done in very small steps to increase accuracy. Once the actual 
simulation starts, the code evaluates for each grid point where the local gas 
is on this temperature evolution and where it should end up for a given 
timestep. It then updates the temperature accordingly so that 

T n A(T re/ ) At 1 



rpn+l _ Y~ 1 



Y{T n ) + 



(16) 



T ref A(T n ) T cooi _ 

where Y is the dimensionless temperature evolution function of a single gas 
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element starting at the maximum temperature of the cooling table as ob- 
tained a priori. In this formula, T re f is an arbitrarily chosen reference tem- 
perature (usually the highest temperature in the cooling table, though this 
isn't requisite). This method has proved to be both faster and more reliable 
than either explicit or implicit methods and is the one that we selected for 
our simulations. For the cool ing curve we use the cooling table calculated by 
Mellema &: Lundqvistl (120021 ) for gas of solar metallicity. The actual (tabu- 
lated) relation between A and T is shown in a log- log plot (fig. [1]), demon- 
strating the extreme challenges due to the large (several orders of magnitude) 
differences. 

The ion density rij and electron density n e follow from the mass density 
and the composition of the gas. For simplicity we assume that hydrogen com- 
pletely dominates the gas (a reasonable assumption for most astrophysical 
situations) and that the gas that actually cools is fully ionized. Therefore, 
Hi = n e = p/rrih, with nth = 1.67 x 10 _24 g the hydrogen mass. 



3.3. Grid setup and AMR strategy 

In order to achieve a high resolution in those areas where it is necessary, 
we use adaptive mesh refinement (AMR). The AMR scheme used in the 
AMRVAC code is a parallellized, mod e rn bl ock-tree variant of the original 
scheme described by iBerger fc Colellal ( 19891 ). This involves the dynamical 
generation (and destruction) of hierarchically nested grids of a fixed subgrid 
(block) size, up to a preset maximum level. This maximum level, along with 
the corresponding effective resolution obtained, is reported for all cases stud- 
ied below. Whether to create or destroy a grid at any giv en level is decided 
automatically at each timestep based on the lLohnerl (119871 ) prescription. This 
estimates a weighted 2 nd derivative of a variable w in grid point i from 



J i'2 i^^^i\^^t2 


d 2 w 


Oxi^ dx^2 


dw 


i-l 


dw 



+ f\w\ 



i+1 



2 ■ 



(17) 



This formula already applies to multi-dimensional simulations, and has a 
wavefilter parameter / which we fix to 10~ 2 . The w indicates an average over 
all neighboring grid cells in directions i x , z 2 . This formula is actually used on 
a user-selected set of variables w from the set of conserved variables U. In our 
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case, we have chosen the density as the variable to initiate refinement. Since 
our simulations are primarily concerned with resolving a high density feature 
(the circumstellar shell), this is the most logical choice. While the formula ITT1 
quantifies the local error for variable w in grid point i, the automated block- 
based regridding works as follows. A fixed gridblock between levels 1 and the 
maximal level allowed is refined if any point on it has this local error exceed a 
pre-set maximum tolerance value set by the user (0.1 in all our simulations). 
In contrast, if all points on the fixed grid-block have their error below a 
fraction (1/8 for our runs) of this maximum tolerance, the grid is coarsened. 
This procedure is then complicated by the proper nesting criterion, ensuring 
no grid changes by more than 1 level across bounding grid blocks. 

3.4- Initial and boundary conditions 

The simulation is set up by using a 1-D spherically symmetric grid with 
400 gridpoints to cover a physical distance starting at 10 17 cm and ending 
at 10 19 cm. Density and pressure in the grid are taken to be constant, with 
the gas velocity set to zero. Boundary conditions at the outer radial bound- 
ary are set to continuous to allow a free outflow. At the inner boundary, 
conditions are set to simulate the inflow of a steady stream of matter ac- 
cording to the values in Table 1. This corresponds to employing standard 
Dirichlet boundary conditions at the inlet. As usual in any finite volume 
code, the boundary conditions include filling of ghost cells, and the stencil of 
our second order shoc k-capturing scheme used demands the use of two ghost 



cell layers. Unlike e.g. Ide Sterck. Rostrup fc Fengl ( 120091 ) we do not attempt 



to simulate the acceleration zone of the wind. Instead, we assume that the 
wind has already reached its terminal velocity (V^) by the time it crosses 
the inner boundary of our grid. The outer boundary is set to continuous (i.e. 
zero gradient extrapolation), though this is by and large irrelevant as the 
simulation is stopped before the expanding wind-bubble reaches the outer 
boundary. 

For our first simulation we use a fixed grid of 400 cells without refinement 
(AMR 1). Subsequently, the simulation is repeated with increasing levels of 
refinement (AMR 2 through AMR 11). Each level of refinement corresponds 
to a doubling of the effective resolution so that AMR 11 corresponds to a 
static grid of 409 600 grid points. To show the effect of radiative cooling on 
both the physical result as well as the numerical simulation we do this for the 
same simulation both with and without radiative cooling. We also include 
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runs with different numerical solvers (van Leer flux-limiter and explicit and 
semi-implicit radiative cooling) to show how they influence the result. 

Finally, we show the results of a single 2-dimensional run to demonstrate 
how shifting to higher dimensions changes the result due to local instabilities 
and to give the reader an idea of what an actual circumstellar shell would 
look like. For this run we use a van Leer flux-limiter, which gives better 
results at higher resolution, with 4 levels of adaptive mesh (as in the 1- 
D AMR 4 simulation). The base grid has 400 radial gridpoints and 200 
azimuthal gridpoints covering the same radius as the 1-D simulations and an 
azimuthal opening angle of 90°. 

3.5. Computational resources 

All calculations were done on the Vic3 supercomputer at K.U.Leuven, 
Belgium, which originally consisted of a 928 core cluster (L5420 CPUs, 1GB 
RAM/core) from SGI, recently updated with 640 cores (Xeon 5560 CPUs, 
3GB RAM/core). All computational nodes are connected through a double 
DDR infiniband network. We typically use a single node of 8 CPUs for 
the 1-D simulations with execution times varying from about 2 minutes for 
the low levels of resolution to a maximum of about 24 hours for the largest 
simulations. 

4. Results 

All simulations described use the setup mentioned earlier. We run simu- 
lations both with and without radiative cooling and compare the results of 
a number of different cooling methods and numerical solvers. We do simula- 
tions with different maximal grid levels (with a factor 2 refinement between 
consecutive levels), compare to uniform grid runs, and vary the physical pa- 
rameters as well. 

4-1. No cooling 

We start with a series of simulations which were run without radiative 
cooling. The results are presented in figs. [2] and [3J Figure |2] shows the 
mass density and temperature of the gas at time 1.25 x 10 12 seconds for 
a simulation with a fixed grid of 400 points. The wind termination shock 
(Rl) and forward shock (R3) are already quite well resolved, although the 
resolution of the forward shock in particular leaves room for improvement. 
The contact discontinuity (R2) is less well resolved and shows the need for a 



13 




_J I I I I I I ! I I I ! I I I I I I ! I I ! LJ -I 

2E+18 4E+18 6E+18 8E+18 

R [cm] 

Figure 2: This figure shows density and temperature after 1.25 x 10 12 seconds for a fixed 
grid simulation with 400 grid points. The wind termination shock (Rl) is well resolved. 
The forward shock (R3) somewhat less so and the contact discontinuity (R2) clearly needs 
a higher resolution. 
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Figure 3: Results for simulations without radiative cooling, showing density (left) and 
temperature (right) after 1.25 x 10 12 seconds (as in fig. [2J. Each simulation has a different 
number of refinement levels, starting from 400 gridpoints with no refinement (AMR1) 
and ending with 4 levels (3 levels on top of the original grid) (AMR4). Note that the 
position of the shell does not change. Only the shape gradually converges, with the 
contact discontinuity R2 most affected by the increased resolution. 



higher effective reso lution. This plot clearly follows the shape predicted by 



Weaver et al.l ( 119771 ). with a free-streaming wind coming from the star, then 
a layer of (near-)constant density shocked wind material (between Rl and 
R2) and finally a swept up shell expanding into t he surrounding inters tellar 
medium. According to the analytical solution by IWeaver et al.l ( 119771 ). the 
outer edge of the shell should be at 6.47 x 10 18 cm for these input parameters. 
In the simulation the outer edge of the shell clearly lies at a larger distance. 
However, the analytical solution depends on the shell being thin compared 
to the shocked wind layer that pushes it outward. In the numerical solution 
this is obviously no t the c ase due to the internal gas-pressure of the shell, 
which IWeaver et al.l ( 119770 neglect. 

Figure [3] shows the same solution as fig. |2] for four simulations with in- 
creasing number of levels of adaptive mesh. The effect of increased resolution 
is small. Both shocks and contact discontinuities become sharper, but the 
shape of the bubble remains the same and the values of density and tem- 
perature are generally the same, except in the transitional regions. The 
simulations with 3 and 4 levels of refinement (AMR 3 and AMR 4 respec- 
tively) overlap each other almost completely, such that a further increase in 
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resolution is unnecessary when pure Eulerian dynamics is studied. 
4-2. With cooling 

The results of our simulations with cooling are given in figs. H] and |5l 
which show the density of the gas at time 1.25 x 10 12 (as in fig. [3]) and time 
2.0 x 10 12 seconds. In fig. H] the simulation without refinement shows clearly 
that the resolution of the grid is too low for an accurate result, as the swept- 
up shell is not properly resolved, looking more like a Gaussian curve than the 
straight angles that should occur at both shocks and discontinuities. Only 
at four levels of refinement and higher the simulations start to give a much 
better result, showing the top of the shell as a flat plateau, though the edges 
are still too rounded and no internal structure in the shell is visible. This is 
solved by adding more resolution. Internal structure of the shell only appears 
at very high levels of resolution (> 6). This structure is a direct result of 
the radiative cooling instability: high density regions cool more rapidly and 
therefore have a lower thermal pressure. As a result they get compressed by 
the surrounding gas, increasing their density even further, which in turn will 
cause them to cool more rapidly. 

Figure [5] shows the same simulations, but at a later time (2.0 x 10 12 sec- 
onds). As time progresses, the amount of mass swept-up in the shell increases, 
which in turn causes the shell to become thicker. As a result it becomes eas- 
ier to resolve. By now even the first simulation is beginning to resolve the 
shell and the 4 level simulation (AMR 4) is already showing internal struc- 
ture. The three highest resolution shells are very similar and overlap each 
other almost completely. Still, the difference in location between the shells at 
low and high refinement shows the need for high effective resolution. When 
comparing figs. H] and for the models with the highest levels of refinement 
we note that the absolute difference does not change over time. The slight 
difference actually originates during the early phases of the simulation. By 
starting the simulation with a direct interaction between wind and stationary 
medium, we created a situation where the initial shocked gas layer will be 
very thin and even a very fine mesh will only be an approximation. Only once 
the shocked gas layer is properly resolved, does the difference in resolution 
no longer matter. 
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Figure 4: Similar to fig. |3l but with radiative cooling. The shell is much thinner due 
to compression. For low resolution the shell is poorly resolved, leading to considerable 
numerical errors. Only for high resolutions, from about AMR 9 onwards, do the results 
converge. We only show the simulations with odd numbered levels of refinement to make 
the figure clearer. 
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Figure 5: Similar to fig. [4j but after 2.0 x f0 12 seconds. The pattern is the same. At 
high resolution the difference between simulations is largely irrelevant, since the absolute 
difference between shell positions stays approximately the same, while the total size of the 
circumstellar bubble increases. 
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Figure 6: The same point in time as in figs. [3] and [4] for two simulations both with 4 levels 
of refinement. One using the exact integration method for the cooling and the other using 
a fully explicit scheme with limited time steps. The results are indistinguishable. 
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Figure 7: Similar to fig. El This time the exact integration method for the cooling is 
compared to a semi-implicit scheme. Again the results are identical. 
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Figure 8: The same point in time as in figs. |3l |4] and [6] for four different simulations: Two 
using the van Leer fiux-limiter method and two using minmod. The van Leer method is 
better at resolving the shell, but at low resolution shows some artificial oscillation inside 
the shocked wind bubble (2 x 10 18 < R < 5 x 10 18 ). 
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Figure 9: The same point in time as in figs. [6] through [8] for two simulations. One with a 
4 level adaptive mesh grid and one with an equal resolution, but a fixed grid. The results 
overlap completely, showing that the adaptive mesh does not introduce numerical errors. 
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Figure 10: The same point in time as in figs. [6] through [9] for 4 simulations with a lower 
mass loss rate and different levels of refinement. The shell has not progressed as far (note 
the different scale on the horizontal axis) and is thicker. Still, the need for refinement 
remains as the internal structure of the shell only becomes visible at high resolution. 
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4-3. Comparing methods 

4-3.1. Different cooling methods 

To make sure that the exact integration method from eq. [THlis equivalent 
to the treatments given by eqs. [T3]and[TI3 we compare the results at the same 
point in time for both the AMR4 simulation with exact integration of the 
cooling curve and identical simulations using fully explicit and semi-implicit 
radiative cooling. For the explicit method the timestep was limited so that 
8e coo i < O.ler, where ct = p/(j — 1) is the internal energy density of the 
gas and 5e coo i is the change in energy density due to radiative cooling. The 
result, shown in fig. El illustrates clearly that there is no noticeable difference 
in the result at all. Similarly, fig. [7] shows the results of the semi-implicit 
versus the exact integration method: again, the results are identical. In this 
particular physical problem, all three methods take approximately equally 
long to compute. The explicit method has the simplest calculation, but oc- 
casionally shorter timesteps are necessary, the exact integration method has 
a longer initial calculation before the actual simulation starts; and the im- 
plicit method requires multiple iterations per timestep. Using one of the 
Vic3 supercomputer's double quad-core Xeon 5560 nodes (2.8Ghz CPUs and 
24GB of RAM) the total computation time for each method is given in Ta- 
ble 2. The main advantage of the exact integration method for this particular 
simulation is its reliability. For the explicit method the maximum allowed 
Se coo i must be chosen in advance and such a choice is by necessity somewhat 
arbitrary. One should actually run multiple simulations to check that the 
timesteps are small enough. The implicit method is vulnerable to potential 
i nstability, due to the particularly complex shape of the cooling curve A(T) 



flTownsendl . l2009f ). 



4-3.2. Flux-limiter influence 

As an alternative to the rather diffusive 'minmod' flux-limiter we also 



tested our cooling routine combined with a 'van Leer' (Ivan Leerl . 119741 1 flux- 
limiter method. Comparing the results for both the AMR 1 and AMR 4 
simulations, an effect becomes apparent (See Fig. [8]). For the lowest reso- 
lution the sharper 'van Leer' limiter yields a result where the shell itself is 
better resolved than with the 'minmod' scheme, but the total bubble is ac- 
tually somewhat smaller. This is due to spurious oscillations in the shocked 
wind bubble (2 x 10 18 < R < 5 x 10 18 ). For the higher resolution, the van 
Leer method yields a result that is clearly much better than the minmod lim- 
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iter. The van Leer AMR 4 result is then reminiscent of a minmod simulation 
with at least one more level of refinement. 

4-3.3. Adaptive mesh versus fixed grid 

To test the effectiveness of the adaptive mesh refinement, we compare 
the AMR 4 result with a fixed grid simulation with a resolution equal to 
the maximum resolution of the adaptive mesh, corresponding to 3200 grid 
points. The result is shown in fig. [HI and appear indistinguishable. From the 
computation time values in Table 2 it is also clear that the adaptive mesh 
scheme saves a considerable amount of time. This is especially important for 
more realistic multi-D simulations. 

4-3.4- Weaker shocks 

In fig. [TU] we show the behaviour of a much weaker shock. For these 
simulations the mass loss rate of the wind has been reduced to 10 -8 M /yr, 
reminiscent of B-type stars (between 2 and 15 M Q ), rather than O-type stars. 
Due to the lower velocity the shock temperatures are much lower (T ~ v 2 ), 
decreasing the effect of the radiative cooling. As a result, the swept-up shell 
which moves at a lower velocity, becomes much more extended. Still, the 
solution does not converge until a relatively high resolution is reached. Al- 
though the wider shell is easier to resolve, the weaker shock will in the earlier 
stages be more radiative (there is less energy available to radiate away, so 
it takes longer to form a shocked wind layer.) This in turn will lead to a 
longer initial phase in which the reaction is almost purely radiative, which is 
far more difficult to resolve. Therefore, although in the later stages the sim- 
ulation will need less resolution, it starts with a larger error. Furthermore, 
although the temperatures in general are lower, the radiative cooling insta- 
bility still occurs, which creates a necessity for a high resolution to resolve 
the internal structure of the shell. 

4-4- Multi-D 

The result of our 2-D simulation is shown in Fig. [TTJ which shows the 
density, temperature and luminosity due to radiative cooling as well as the 
refinement of the grid. For this simulation we use a van Leer flux-limiter, 
which gives a better performance at lower resolution (see fig. ED, and we 
allow a maximum of 4 levels of refinement. In two dimensions the shell is 
not a perfect spherical shape. Rather, it shows considerable structure. This 
is primarily due to Ray leigh- Taylor instability, induced by the fact that a 



25 



iog(L) 



8E+18 - 



4E+18 - 



-24 -22.25 -20.5 

1 I 1 1 1 1 I 




-8E+18 - 



-8E+18 



8E+18 



Figure 11: The result of a 2-dimensional simulation, showing (starting upper left and 
moving clockwise) the luminosity, density, temperature and grid-structure. All units are 
in cgs format. The structure of the shell is quite complicated due to Ray leigh- Taylor 
instability. The adaptive mesh refinement has increased the grid resolution around the 
areas with steep density jumps (The shell and the reverse shock). The luminosity, which 
is strongly dependent on the density, only shows the shell itself. 



26 



Table 2: Computation time in CPU seconds 

method total time (s) Startup time (s) I/O (s 

exact integration 232.925 0.081 0.394 

explicit 243.365 0.013 0.444 

semi-implicit 253.205 0.038 0.699 

exact integration (3200 points, 509.003 0.072 0.829 
no AMR) 



low density bubble (2 x 10 18 < R < 5 x 10 18 ) expands into a much denser 
medium. The luminosity plot only shows the shell, owing to the strong 
density dependence of the radiative cooling (see eq. [3]). The high temperature 
bubble of shocked wind material barely radiates on account of its low density. 
The adaptive mesh has refined mostly around the shell and the reverse shock. 
We stress that the plot does not show the actual gridcells, these would be 
too small to see. What is shown is the fixed size subgrids, each consisting of 
20x20 gridcells. 

5. Discussion of astrophysical relevance 

Our results have shown the increased need for high resolution caused by 
the addition of radiative cooling to a numerical simulation. Since astrophysics 
usually deals with phenomena on very large scales, adaptive mesh refinement 
becomes a prerequisite in order to keep the total size of the simulation to 
a manageable level. Of course, this particular test-case does not cover all 
possible physical processes that contribute to an astrophysical problem. E.g. 
we have neglected the issue of an existing radiation field that may heat the 
local gas. Nor did we include thermal conduction, or magnetic fields. In 
reality, all these contribute, but it falls outside the scope of this paper to 
address them all. Nevertheless, we can make qualitative predictions of their 
influence on the final result. 

5.1. Magnetic fields 

Magnetic fields play an important role in the evolution of the circumstel- 
lar medium. Not only do magnetic fields exist in the interstellar gas, but the 
stars themselves can have powerfull magnetic fields that influence the char- 
acteristics of the stellar wind. The magnetic field strength in the interstellar 
medium is known to vary from ~ 1/zG in the general interstellar medium 
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flBeckl . 120091 : Ijaffe et all l2010h to 0.1-6 mG (ICurran fc Chrvsostomoul l2007t ) 
in some massive star forming regions. As such, the magnetic field energy 
density is generally low, although this can increase if the gas is swept up by 
the stellar wind. Much will depend on the alignment of the magnetic field 
lines with the motion of the swept up shell. If the field lines run parallel with 
the direction in which the shell moves, the influence of the field on the general 
morphology of the gas is likely negligible. However, should the field lines run 
parallel with the shell, the local magnetic field energy density increases as 
the gas is compressed. This causes the shell to expand as the magnetic field 
counteracts the pressure, resulting in a lower density of the shell. This acts 
to decrease the effect of radiative cooling. 

The magnetic field of the star itself is unlikely to be rele yant, although 



some massiv e stars have strong magnetic fields. Simulations by lud-Doula. Owocki & Townsend 
(120081 . 120091 ) show that far from the star (> 10 stellar radii) the magnetic 
field lines are torn open by the wind. This results in a monopolar field struc- 
ture, where the field lines become parallel with the direction of the wind. 
Therefore, the stellar field will not influence the kinetics of the circumstellar 
gas at large distances. 

5.2. Radiation fields 

The absorption of high energy photons from an outside source (such as 
nearby stars) can influence the morphology of the circumstellar medium and 
photo-ionize the circumstellar g as. The effect of photo-ionization was de 



scrib ed in a simple ID model by Ivan Marie. Langer &: Garcia-Segural (12 005 



2007) a nd 2D simulations of similar problems were run by lFreyer. Hensler. fc Yorke 
( 120031 . 120061 ) . The main effect is the increased temperature of the ma- 
terial into which the stellar winds expands. This raises the sound speed 
and depending on the strength of the wind, the expansion speed of the 
swept-up shell can drop below this sound speed. Then, the forward shock 
is lost and the transition between expanding shocked wind and surround- 
i ng (ionized) medium becomes a contact d iscontinuity rather than a shock 
( Ivan Marie. Langer fc Garcia-Segural . 120051 ) . At the outer edge of the photo- 
ionized region, the thermal pressure of the photo-ionized gas will still cause a 
shock to form as the ionized gas expands into the cold surrounding medium. 
On the other hand, if the wind is strong enough to maintain a supersonic 
expansion, the result will be almost identical to the si tuation without photo- 
ionization (Ivan Marie. Langer fe Garcia-Segural . 120071 ). Photo-ionization can 
also increase the instability of a partially ionized shell, as local photo-ionized 



28 



region s have a much higher t empe rature than their surroundings and will ex- 
pand ( Garcia-Segura. et al. . 19991 ) . However, the shells we find already have 



a temperature of about 10,000 K (the ionization temperature of hydrogen) 
due to collisional heating. Therefore, it is unlikely that ionization due to the 
presence of radiation will have a large effect on the structure of the shell. 

5.3. Going to 3-D 

As shown in section I4.4[ the circumstellar shell is highly susceptible to 
instabilities in more than 1 spatial dimension. The full implications of thin- 
shell instabilities in 3D have not yet been explored in great detail, and form 
part of our ongoing efforts. Generally though, instabilities tend to form 



quicker in 3D than in 2D (e.g. lYoung et al.l . 120011 ). quickly producing higher 



density contrasts and creating more complicated structures. This can only 
increase the effects described in this paper, as higher densities lead inevitably 
to greater radiative cooling, which in turn will increase the density as the 
cooling gas is compressed. Once again, there is a strong need for a high 
resolution grid and therefore the application of AMR will be inevitable. 

6. Conclusion 

In order to properly resolve the thin, high density shells resulting from ra- 
diative cooling, astrophysical gas dynamics simulations need grids with high 
resolution. Since the number of gridpoints would quickly become prohibitive 
if this simulation was carried out on a fixed grid, adaptive mesh refinement 
can be considered a necessity for simulations of this kind. In 2-D (and 3-D) 
the adaptive mesh becomes even more crucial, since the number of gridpoints 
is much larger to begin with and the complicated structures need high reso- 
lution to be properly resolved. Numerically, a better result could be achieved 
by replacing the diffusive minmod scheme in a linear extrapolation with a 
higher order interpo l ation method, such as the piec ewise parabolic method 



(IFrvxell et all 1200(1 Mignone. Plewa fc Bodol . 120051 1 . This would allow the 



code to reach a more accurate solution at lower resolution. 

Finally, this particular simulation has rather high densities and shock 
temperature due to the high mass loss rate of the stellar wind. This in 
turns increases the effect of the radiative cooling, due to its strong density 
dependence. Weaker shocks can be equally problematic to resolve due to their 
more radiative nature in the inital collision. This could perhaps be mitigated 
by starting the simulation gradually, rather than have the full strength wind 
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interact directly with the stationary gas. Unfortunately, any simulation then 
depends on parameters used to control the transition. Also, high density 
regions will still be compressed due to radiative cooling, which means that 
high resolution grids will be needed. 

We also stress the actual scale that these simulations cover. 409 600 grid- 
cells along one dimension may seem excessive, but since the whole grid covers 
a range of 10 19 cm, this means that each gridcell has a cross-section of about 
35 solar radii. Of course, one can wonder whether a proper resolution of the 
shell is really necessary. For one thing, typical astrophysical problems involve 
significant uncertainties in the input parameters, which dwarf the possible 
effect of the resolution. However, for anyone working on such simulations 
it is prerequisite to be aware of this particular problem and consider how 
it may affect scientific results. A viable alternative to AMR grids as used 
here, is the u se of moving grids, whi c h hay e recently been used in pure 2D 
Euler flows in Ivan Dam and Zegelingl (j2Q10h . The effect of radiative losses as 
emphasized here still need to be evaluated for moving grid simulations. 
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